Create datasets

# Load data on individual donations
load("./data/r02.fd.bd.all.rdata") # outputs an object called "output" into the environment
donations <- output

# We only want to look at first donation event values from each donor
donors <- donations %>%
    group_by(donor) %>%
    filter(date == min(date)) %>%
    ungroup()

# Load FinDonor demographic data
load("./data/r02ds.donorData.rdata") # outputs an object called "output" into the environment
findonor <- output

# Combine the FinDonor datasets
FinDonor <- left_join(donors, findonor, by = "donor")

# Load THL data
# Sofie: thldalta.rdata contains all five THL cohorts, extract FINRISK97 and Health2000 from the others
load("./data/thldata.rdata")
FinRisk97 <- thldata$fr1997
Health2k <- thldata$h2000

# Remove leftovers
rm(output)
rm(thldata)

## Rename useful stuff
# Ferritin, Self-Reported Health
FinDonor <- rename(FinDonor, SRH = QR17, Menstruation = QR79, Age_float = Age, Age = age)
FinRisk97 <- rename(FinRisk97, Ferritin = FERRITIN, SRH = Q40, Gender = SUKUP, Menstruation = K129, Age = IKA)
Health2k <- rename(Health2k, Ferritin = FERRITIINI, SRH = BA01, Gender = SP2, Menstruation = BD03, Age = IKA2, Menopause = MENOP, APOB = NMR_APOB, APOA1 = NMR_APOA1)

# Make "useful stuff" conform with each other
FinDonor <- FinDonor %>% 
    mutate(SRH = case_when(SRH == "Excellent" ~ 1,
                           SRH == "Very_good" ~ 2,
                           SRH == "Good" ~ 3,
                           SRH == "Satisfactory" ~ 4,
                           SRH == "Poor" ~ 5),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == "regular_period" | Menstruation == "irregular_period") ~ "Women|Menstr",
                             Gender == "Women" & Menstruation == "no_period" ~ "Women|Non-menstr",
                             TRUE ~ "NA")) # Equates to "else"

FinRisk97 <- FinRisk97 %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Menstr",
                             Gender == "Women" & Menstruation == 3 ~ "Women|Non-menstr",
                             TRUE ~ "NA"))

Health2k <- Health2k %>%
    mutate(Gender = case_when(Gender == 1 ~ "Men",
                              Gender == 2 ~ "Women",
                              TRUE ~ "NA"),
           # It is worth noting here, that menstruation status was not asked from women over 55. So, we treat these as postmenopausal.
           Group = case_when(Gender == "Men" ~ "Men",
                             Gender == "Women" & (Menstruation == 1 | Menstruation == 2) ~ "Women|Menstr",
                             Gender == "Women" & (Menstruation == 3 | Age >= 55) ~ "Women|Non-menstr",
                             TRUE ~ "NA"))

# Donation eligibility
# These are both "approximates" in a sense, we don't have all the necessary variables to
# filter thoroughly, and we'll be able to do more filtering on Health2000 than FinRisk97
donor_eligible_h2k <- Health2k %>% 
    filter(BMII_PAINO.x >= 50 & BMII_PAINO.x <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 & Age <= 66) %>% # Filter away too young and too old
    filter((B_Hb >= 125 & Gender == "Women") | (B_Hb >= 135 & Gender == "Men")) %>% # Filter by hemoglobin
    filter(BA08 == 0) %>% # filter out people with heart attacks
    filter(BA09 == 0) %>% # filter out people with angina
    filter(BA10 == 0) %>% # cardiac insufficiency / heart failure
    filter(!(BA26 == 1 & ATC_A10A == 1)) %>% # filter out people who are diabetic AND use insulin
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) %>% # rename for ease of use
    mutate(HbA1C = B_GHb_A1C * 10.93 - 23.50)

donor_eligible_fr <- FinRisk97 %>%
    filter(PAINO >= 50 & PAINO <= 200) %>% # Filter away people <50kg and >200kg
    filter(Age >= 18 & Age <= 66) %>% # Filter away too young and too old
    #filter((HGB >= 125 & Gender == 2) | (HGB >= 135 & Gender == 1)) %>% # DON'T filter by hemoglobin, < 500 values in data
    filter(Q15A != 2) %>% # STEMI, NSTEMI
    filter(Q16A != 2) %>% # Stroke
    # filter(Q38 != 2 & Q38 != 4) %>%  # Insulin treatment (2: just insulin, 4: insulin and a tablet)
    filter(Q17B != 2) %>% # cardiac insufficiency
    filter(Q17C != 2) %>% # angina pectoris
    filter(SRH < 4) %>% # filter out "Bad" or "Very bad" SRH
    rename(GlycA = GP) # rename for ease of use

# Create useful mastersets
# fer_srh <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "SRH", "Group")], 
#                      Health2k = donor_eligible_h2k[, c("Ferritin", "SRH", "Group")], .id = "Cohort") %>% 
#     mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
#            Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k")),
#            SRH = ordered(SRH, levels = 1:5)) %>%
#     filter(Group != "NA") %>%
#     drop_na()

fer_crp <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "CRP")], 
                         Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "CRP")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    filter(CRP >= 0.01) %>%
    drop_na()

Mock Cohort|Var Tables

table1 <- as.data.frame(table(fer_crp$Group, fer_crp$Cohort))
table1$CRP <- c(paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$CRP[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1$FER <- c(paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Menstr" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Women|Non-menstr" & fer_crp$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_crp$Ferritin[fer_crp$Group == "Men" & fer_crp$Cohort == "Health2k"])[5], 2), ")"))
table1
##               Var1      Var2 Freq                 CRP                      FER
## 1     Women|Menstr FinRisk97 1915  0.77 | (0.4, 1.87)   24.02 | (12.39, 42.53)
## 2 Women|Non-menstr FinRisk97  884 1.27 | (0.62, 2.65)   55.79 | (31.06, 93.72)
## 3              Men FinRisk97 2603 0.89 | (0.46, 1.88) 112.06 | (66.46, 182.07)
## 4     Women|Menstr  Health2k  943 0.62 | (0.27, 1.81)        28 | (15.2, 48.6)
## 5 Women|Non-menstr  Health2k  661  1.02 | (0.38, 2.4)      55.9 | (32.7, 95.5)
## 6              Men  Health2k 1710 0.77 | (0.35, 1.75)   124.7 | (76.67, 193.9)

Exclusions: Health2000

# ExclusionTable method
excl_h2k_women_mens <- exclusion_table(data = Health2k,
                inclusion_criteria = c("Group == 'Women|Menstr'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "BMII_PAINO.x < 50", "BMII_PAINO.x > 200", "Age < 18",
                                       "Age > 66", "SRH >= 4", "B_Hb < 125", "BA08 == 1", "BA09 == 1", "BA10 == 1", "BA26 == 1 & ATC_A10A == 1"),
                keep_data = TRUE)
excl_h2k_women_mens
## 
## ======================================================
## Excluded the following observations:
## ======================================================
## Exclusions based on INCLUSION criteria
## 
##                 inclusion n_prior n_post n_excluded
## 1 Group == 'Women|Menstr'    6264   1508       4756
## 2                   TOTAL    6264   1508       4756
## 
## Exclusions based on EXCLUSION criteria
## 
##                    exclusion n_prior n_post n_excluded
## 1                 is.na(CRP)    1508   1464         44
## 2            is.na(Ferritin)    1464   1442         22
## 3                 CRP < 0.01    1442   1270        172
## 4          BMII_PAINO.x < 50    1270   1230         40
## 5         BMII_PAINO.x > 200    1230   1230          0
## 6                   Age < 18    1230   1230          0
## 7                   Age > 66    1230   1230          0
## 8                   SRH >= 4    1230   1187         43
## 9                 B_Hb < 125    1187    953        234
## 10                 BA08 == 1     953    952          1
## 11                 BA09 == 1     952    951          1
## 12                 BA10 == 1     951    949          2
## 13 BA26 == 1 & ATC_A10A == 1     949    943          6
## 14                     TOTAL    1508    943        565
## 
## ======================================================
excl_h2k_women_nonmens <- exclusion_table(data = Health2k,
                inclusion_criteria = c("Group == 'Women|Non-menstr'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "BMII_PAINO.x < 50", "BMII_PAINO.x > 200", "Age < 18",
                                       "Age > 66", "SRH >= 4", "B_Hb < 125", "BA08 == 1", "BA09 == 1", "BA10 == 1", "BA26 == 1 & ATC_A10A == 1"),
                keep_data = TRUE)
excl_h2k_women_nonmens
## 
## ========================================================
## Excluded the following observations:
## ========================================================
## Exclusions based on INCLUSION criteria
## 
##                     inclusion n_prior n_post n_excluded
## 1 Group == 'Women|Non-menstr'    6264   1505       4759
## 2                       TOTAL    6264   1505       4759
## 
## Exclusions based on EXCLUSION criteria
## 
##                    exclusion n_prior n_post n_excluded
## 1                 is.na(CRP)    1505   1160        345
## 2            is.na(Ferritin)    1160   1148         12
## 3                 CRP < 0.01    1148   1078         70
## 4          BMII_PAINO.x < 50    1078   1061         17
## 5         BMII_PAINO.x > 200    1061   1061          0
## 6                   Age < 18    1061   1061          0
## 7                   Age > 66    1061    878        183
## 8                   SRH >= 4     878    782         96
## 9                 B_Hb < 125     782    706         76
## 10                 BA08 == 1     706    700          6
## 11                 BA09 == 1     700    678         22
## 12                 BA10 == 1     678    668         10
## 13 BA26 == 1 & ATC_A10A == 1     668    661          7
## 14                     TOTAL    1505    661        844
## 
## ========================================================
excl_h2k_men <- exclusion_table(data = Health2k,
                inclusion_criteria = c("Group == 'Men'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "BMII_PAINO.x < 50", "BMII_PAINO.x > 200", "Age < 18",
                                       "Age > 66", "SRH >= 4", "B_Hb < 135", "BA08 == 1", "BA09 == 1", "BA10 == 1", "BA26 == 1 & ATC_A10A == 1"),
                keep_data = TRUE)
excl_h2k_men
## 
## ======================================================
## Excluded the following observations:
## ======================================================
## Exclusions based on INCLUSION criteria
## 
##        inclusion n_prior n_post n_excluded
## 1 Group == 'Men'    6264   2944       3320
## 2          TOTAL    6264   2944       3320
## 
## Exclusions based on EXCLUSION criteria
## 
##                    exclusion n_prior n_post n_excluded
## 1                 is.na(CRP)    2944   2517        427
## 2            is.na(Ferritin)    2517   2494         23
## 3                 CRP < 0.01    2494   2254        240
## 4          BMII_PAINO.x < 50    2254   2251          3
## 5         BMII_PAINO.x > 200    2251   2251          0
## 6                   Age < 18    2251   2251          0
## 7                   Age > 66    2251   2103        148
## 8                   SRH >= 4    2103   1886        217
## 9                 B_Hb < 135    1886   1791         95
## 10                 BA08 == 1    1791   1756         35
## 11                 BA09 == 1    1756   1735         21
## 12                 BA10 == 1    1735   1727          8
## 13 BA26 == 1 & ATC_A10A == 1    1727   1710         17
## 14                     TOTAL    2944   1710       1234
## 
## ======================================================

Exclusions: FinRisk1997

# ExclusionTable method
excl_fr_women_mens <- exclusion_table(data = FinRisk97,
                inclusion_criteria = c("Group == 'Women|Menstr'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "PAINO < 50", "PAINO > 200", "Age < 18", 
                                       "Age > 66", "SRH >= 4", "Q15A == 2", "Q16A == 2", "Q17B == 2", "Q17C == 2"),
                keep_data = TRUE)
excl_fr_women_mens
## 
## ====================================================
## Excluded the following observations:
## ====================================================
## Exclusions based on INCLUSION criteria
## 
##                 inclusion n_prior n_post n_excluded
## 1 Group == 'Women|Menstr'    7943   2469       5474
## 2                   TOTAL    7943   2469       5474
## 
## Exclusions based on EXCLUSION criteria
## 
##          exclusion n_prior n_post n_excluded
## 1       is.na(CRP)    2469   2288        181
## 2  is.na(Ferritin)    2288   2159        129
## 3       CRP < 0.01    2159   2159          0
## 4       PAINO < 50    2159   2092         67
## 5      PAINO > 200    2092   2092          0
## 6         Age < 18    2092   2092          0
## 7         Age > 66    2092   2084          8
## 8         SRH >= 4    2084   1977        107
## 9        Q15A == 2    1977   1971          6
## 10       Q16A == 2    1971   1962          9
## 11       Q17B == 2    1962   1939         23
## 12       Q17C == 2    1939   1915         24
## 13           TOTAL    2469   1915        554
## 
## ====================================================
excl_fr_women_nonmens <- exclusion_table(data = FinRisk97,
                inclusion_criteria = c("Group == 'Women|Non-menstr'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "PAINO < 50", "PAINO > 200", "Age < 18", 
                                       "Age > 66", "SRH >= 4", "Q15A == 2", "Q16A == 2", "Q17B == 2", "Q17C == 2"),
                keep_data = TRUE)
excl_fr_women_nonmens
## 
## ========================================================
## Excluded the following observations:
## ========================================================
## Exclusions based on INCLUSION criteria
## 
##                     inclusion n_prior n_post n_excluded
## 1 Group == 'Women|Non-menstr'    7943   1463       6480
## 2                       TOTAL    7943   1463       6480
## 
## Exclusions based on EXCLUSION criteria
## 
##          exclusion n_prior n_post n_excluded
## 1       is.na(CRP)    1463   1389         74
## 2  is.na(Ferritin)    1389   1319         70
## 3       CRP < 0.01    1319   1319          0
## 4       PAINO < 50    1319   1289         30
## 5      PAINO > 200    1289   1289          0
## 6         Age < 18    1289   1289          0
## 7         Age > 66    1289   1145        144
## 8         SRH >= 4    1145   1002        143
## 9        Q15A == 2    1002    986         16
## 10       Q16A == 2     986    959         27
## 11       Q17B == 2     959    908         51
## 12       Q17C == 2     908    884         24
## 13           TOTAL    1463    884        579
## 
## ========================================================
excl_fr_men <- exclusion_table(data = FinRisk97,
                inclusion_criteria = c("Group == 'Men'"),
                exclusion_criteria = c("is.na(CRP)", "is.na(Ferritin)", "CRP < 0.01", "PAINO < 50", "PAINO > 200", "Age < 18", 
                                       "Age > 66", "SRH >= 4", "Q15A == 2", "Q16A == 2", "Q17B == 2", "Q17C == 2"),
                keep_data = TRUE)
excl_fr_men
## 
## ============================================
## Excluded the following observations:
## ============================================
## Exclusions based on INCLUSION criteria
## 
##        inclusion n_prior n_post n_excluded
## 1 Group == 'Men'    7943   3943       4000
## 2          TOTAL    7943   3943       4000
## 
## Exclusions based on EXCLUSION criteria
## 
##          exclusion n_prior n_post n_excluded
## 1       is.na(CRP)    3943   3685        258
## 2  is.na(Ferritin)    3685   3476        209
## 3       CRP < 0.01    3476   3476          0
## 4       PAINO < 50    3476   3469          7
## 5      PAINO > 200    3469   3469          0
## 6         Age < 18    3469   3469          0
## 7         Age > 66    3469   3159        310
## 8         SRH >= 4    3159   2833        326
## 9        Q15A == 2    2833   2773         60
## 10       Q16A == 2    2773   2729         44
## 11       Q17B == 2    2729   2652         77
## 12       Q17C == 2    2652   2603         49
## 13           TOTAL    3943   2603       1340
## 
## ============================================

Ferritin X CRP

Subgroups: menstruating women, non-menstruating women, men. Grey box highlights all individuals that go over 3 mg/l CRP at ferritin >15 ug/l.

options(scipen = 10000)
ratio1 <- round(nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97")), 2)
ratio3 <- round(nrow(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97")), 2)
ratio4 <- round(nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 6), CRP = rep(50, 6), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%"), paste0(ratio3, "%"), paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = c(rep("FinRisk97", 3), rep("Health2k", 3)), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 3, ymax = 200, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    theme_minimal() + 
    geom_smooth(method = "lm", formula = y ~ x, color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "CRP (mg/l)")

Ferritin X CRP | COLORED

options(scipen = 10000)
ratio1 <- round(nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97")), 2)
ratio3 <- round(nrow(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97")), 2)
ratio4 <- round(nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 6), CRP = rep(50, 6), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%"), paste0(ratio3, "%"), paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = c(rep("FinRisk97", 3), rep("Health2k", 3)), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 3, ymax = 200, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none") +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "CRP (mg/l)")

Ferritin X CRP | GROUPS ONLY BY COLOR

options(scipen = 10000)
ratio1 <- round(nrow(fer_crp %>% filter(Cohort == "FinRisk97", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_crp %>% filter(Cohort == "Health2k", CRP >= 3 & Ferritin >= 15)) * 100 / nrow(fer_crp %>% filter(Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 2), CRP = rep(50, 2), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%")), Cohort = c("FinRisk97", "Health2k"))
ggplot(data = fer_crp, aes(x = Ferritin, y = CRP)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 3, ymax = 200, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.2) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(cols = vars(Cohort)) +
    theme(legend.position = "bottom") + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "CRP (mg/l)")

Proportion analysis code

The estimate: How many more people are over 3 mg/l in CRP when we filter the populations by different ferritin thresholds? Expressed in percentage point difference. The intuition: We currently don’t do ferritin filtering. The population from which we draw donors has a certain proportion of people in risk of cardiac/chronic inflammation, often indicated by elevated hs-CRP (>3 mg/l). Does this proportion increase in a statistically significant manner, if we start requiring higher ferritin levels from donors?

if (!file.exists(paste0("./data/PUBL_finrisk_CRP_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        bootobs$men[[i]] <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        # Women|Menstr
        bootobs$women_mens[[i]] <- boot(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        # Women|Non-menstr
        bootobs$women_nonmens[[i]] <- boot(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i],var2_trld = CRP_trld, var2_over = T)
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        
        ## Store
        # Men
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Menstr
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        # Women|Non-menstr
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_mens, means_women_nonmens),
                                     upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                     lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                     Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_finrisk_CRP_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_finrisk_CRP_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    
    # Combine
    means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

if (!file.exists(paste0("./data/PUBL_health2k_CRP_", boot_n, ".rds"))) { # run bootstrap only if needed
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    bootobs$men[[i]] <- boot(fer_crp %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
    # Women|Menstr
    bootobs$women_mens[[i]] <- boot(fer_crp %>% filter(Group == "Women|Menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
    # Women|Non-menstr
    bootobs$women_nonmens[[i]] <- boot(fer_crp %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = CRP, var1_trld = ferritin_values[i], var2_trld = CRP_trld, var2_over = T)
    ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
    
    ## Store
    # Men
    means_men[i] <- bootobs$men[[i]]$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Menstr
    means_women_mens[i] <- bootobs$women_mens[[i]]$t0
    upper_women_mens[i] <- ci_obj_women_mens$normal[3]
    lower_women_mens[i] <- ci_obj_women_mens$normal[2]
    # Women|Non-menstr
    means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
    upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
    lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))    
    
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_health2k_CRP_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_health2k_CRP_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 153), rep("Health2k", 153))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Menstr", "Women|Non-menstr", "Men"))

Proportion analysis plot

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") + guides(linetype = "none")

Checking statistical significance at points of interest

IMPORTANT: We are checking for non-overlaps of 95% confidence intervals. These intervals have been derived from 10000 bootstrap samples using normal approximation. Normal approximation relies on the assumption, that our bootstrapped estimates are normally distributed, which we have confirmed separately using both visual assessment (a histogram) and a standard Shapiro-Wilk test.

# FINRISK
# Menstruating Women
frwomenpre5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Menstr" & Ferritin == 5)
frwomenpre15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Menstr" & Ferritin == 15)
frwomenpre30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Menstr" & Ferritin == 30)
frwomenpre50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Menstr" & Ferritin == 50)

# Non-menstruating Women
frwomenpost5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Non-menstr" & Ferritin == 5)
frwomenpost15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Non-menstr" & Ferritin == 15)
frwomenpost30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Non-menstr" & Ferritin == 30)
frwomenpost50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Women|Non-menstr" & Ferritin == 50)

# Men
frmen5 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 5)
frmen15 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 15)
frmen30 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 30)
frmen50 <- means_all %>% filter(Cohort == "FinRisk97" & Group == "Men" & Ferritin == 50)

# H2K
# Menstruating Women
h2kwomenpre5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Menstr" & Ferritin == 5)
h2kwomenpre15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Menstr" & Ferritin == 15)
h2kwomenpre30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Menstr" & Ferritin == 30)
h2kwomenpre50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Menstr" & Ferritin == 50)

# Non-menstruating Women
h2kwomenpost5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Non-menstr" & Ferritin == 5)
h2kwomenpost15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Non-menstr" & Ferritin == 15)
h2kwomenpost30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Non-menstr" & Ferritin == 30)
h2kwomenpost50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Women|Non-menstr" & Ferritin == 50)

# Men
h2kmen5 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 5)
h2kmen15 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 15)
h2kmen30 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 30)
h2kmen50 <- means_all %>% filter(Cohort == "Health2k" & Group == "Men" & Ferritin == 50)

Using these objects, we can check the significance of between the points of interest. For example “h2kwomenpre5\(upper < h2kwomenpre15\)lower” would evaluate to TRUE, if the proportion of menstruating women over 3 mg/l CRP is significantly higher in the population filtered by ferritin 15 ug/l than 5 ug/l. Indeed, for menstruating women, the proportion is significantly higher at filter levels of 15, 30, and 50 ug/l when compared with 5 ug/l. This holds for non-menstruating women also, except for the FinRisk 1997 cohort, where the difference between levels 5 and 15 was not significant. The differences are all significant also in men, but the respective increases in proportions are much smaller.

Proportion analysis plot | COLORED

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") +
    theme(legend.position = "none")

Proportion analysis plot | COLORED & ALL-IN-1

ggplot(data = means_all, aes(x = Ferritin, y = means, group = Group)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    #geom_point(aes(color = Group, shape = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") +
    guides(fill = "none") +
    theme(legend.position = "bottom")

Additional markers (for the supplement)

We’ll also take a look at acetylated glycoprotein measurements, glucated haemoglobin, and apolipoproteins A1 and B.

# mastersets for the Supplement
# GlycA
fer_glyca <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "GlycA")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "GlycA")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# HbA1C
fer_hba1c <- bind_rows(Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "HbA1C")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
           Cohort = ordered(Cohort, levels = c("Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()

# APOB
fer_apob <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOB")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOB")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()
# APOA1
fer_apoa1 <- bind_rows(FinRisk97 = donor_eligible_fr[, c("Ferritin", "Group", "APOA1")], 
                       Health2k = donor_eligible_h2k[, c("Ferritin", "Group", "APOA1")], .id = "Cohort") %>% 
    mutate(Group = ordered(Group, levels = c("Women|Menstr", "Women|Non-menstr", "Men")),
           Cohort = ordered(Cohort, levels = c("FinRisk97", "Health2k"))) %>%
    filter(Group != "NA") %>%
    drop_na()

Mock Cohort|Var Tables (Supplement)

suptable1 <- as.data.frame(table(fer_glyca$Group, fer_glyca$Cohort))

suptable1$GlycA <- c(paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$GlycA[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable1$FER <- c(paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Menstr" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Women|Non-menstr" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_glyca$Ferritin[fer_glyca$Group == "Men" & fer_glyca$Cohort == "Health2k"])[5], 2), ")"))

suptable2 <- as.data.frame(table(fer_hba1c$Group))
suptable2$HbA1C <- c(paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$HbA1C[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable2$FER <-  c(paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Menstr" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Women|Non-menstr" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_hba1c$Ferritin[fer_hba1c$Group == "Men" & fer_hba1c$Cohort == "Health2k"])[5], 2), ")"))

suptable1
##               Var1      Var2 Freq               GlycA                     FER
## 1     Women|Menstr FinRisk97 1977 1.28 | (1.17, 1.42)    23.9 | (12.3, 42.52)
## 2 Women|Non-menstr FinRisk97  891 1.38 | (1.25, 1.52)   55.78 | (31.2, 93.21)
## 3              Men FinRisk97 2640 1.38 | (1.25, 1.55)  112.1 | (66.53, 183.4)
## 4     Women|Menstr  Health2k 1072  1.1 | (0.96, 1.23)  27.64 | (14.78, 48.25)
## 5 Women|Non-menstr  Health2k  704     1.15 | (1, 1.3)      55.6 | (32, 95.25)
## 6              Men  Health2k 1903 1.19 | (1.05, 1.35) 121.8 | (75.72, 189.34)
suptable2
##               Var1 Freq                  HbA1C                    FER
## 1     Women|Menstr 1074 31.15 | (28.96, 33.34) 27.62 | (14.72, 48.19)
## 2 Women|Non-menstr  705 33.34 | (31.15, 36.61)      55.6 | (32, 95.4)
## 3              Men 1910 34.43 | (32.24, 35.52) 121.8 | (75.8, 189.42)
suptable3 <- as.data.frame(table(fer_apob$Group, fer_apob$Cohort))

suptable3$APOB <- c(paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$APOB[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[5], 2), ")"))

suptable3$FER <- c(paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Menstr" & fer_apob$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Women|Non-menstr" & fer_apob$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apob$Ferritin[fer_apob$Group == "Men" & fer_apob$Cohort == "Health2k"])[5], 2), ")"))

suptable3
##               Var1      Var2 Freq                APOB                     FER
## 1     Women|Menstr FinRisk97 1935 0.86 | (0.73, 1.01)  24.05 | (12.38, 42.51)
## 2 Women|Non-menstr FinRisk97  886   1.04 | (0.9, 1.2)  55.79 | (31.14, 94.17)
## 3              Men FinRisk97 2615  1.01 | (0.85, 1.2) 111.93 | (65.81, 181.7)
## 4     Women|Menstr  Health2k 1072 0.75 | (0.63, 0.89)  27.64 | (14.78, 48.25)
## 5 Women|Non-menstr  Health2k  704 0.88 | (0.74, 1.03)      55.6 | (32, 95.25)
## 6              Men  Health2k 1904  0.93 | (0.78, 1.1) 121.8 | (75.72, 189.57)
suptable4 <- as.data.frame(table(fer_apoa1$Group, fer_apoa1$Cohort))

suptable4$APOA1 <- c(paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$APOA1[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"))

suptable4$FER <- c(paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "FinRisk97"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Menstr" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Women|Non-menstr" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"),
                  paste0(round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[3], 2), " | (", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[2], 2), ", ", round(summary(fer_apoa1$Ferritin[fer_apoa1$Group == "Men" & fer_apoa1$Cohort == "Health2k"])[5], 2), ")"))

suptable4
##               Var1      Var2 Freq               APOA1                      FER
## 1     Women|Menstr FinRisk97 1935   1.7 | (1.51, 1.9)   24.05 | (12.38, 42.51)
## 2 Women|Non-menstr FinRisk97  886 1.73 | (1.53, 1.93)   55.79 | (31.14, 94.17)
## 3              Men FinRisk97 2614  1.5 | (1.34, 1.68) 111.72 | (65.77, 181.74)
## 4     Women|Menstr  Health2k 1072 1.51 | (1.35, 1.66)   27.64 | (14.78, 48.25)
## 5 Women|Non-menstr  Health2k  704 1.55 | (1.38, 1.72)       55.6 | (32, 95.25)
## 6              Men  Health2k 1904 1.44 | (1.31, 1.57)  121.8 | (75.72, 189.57)

Ferritin X GlycA | COLORED

options(scipen = 10000)
ratio1 <- round(nrow(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97")), 2)
ratio3 <- round(nrow(fer_glyca %>% filter(Group == "Men" & Cohort == "FinRisk97", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Men" & Cohort == "FinRisk97")), 2)
ratio4 <- round(nrow(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_glyca %>% filter(Group == "Men" & Cohort == "Health2k", GlycA >= 1.35 & Ferritin >= 15)) * 100 / nrow(fer_glyca %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 6), GlycA = rep(2.2, 6), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%"), paste0(ratio3, "%"), paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = c(rep("FinRisk97", 3), rep("Health2k", 3)), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_glyca, aes(x = Ferritin, y = GlycA)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 1.35, ymax = 3.7, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none") +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")),
         y = "GlycA (mmol/l)")

Proportion analysis code (GlycA)

For the acetylated glycoprotein we don’t have any well established thresholds for “healthy” and “unhealthy”. We’ll use the population median here.

iterations <- length(ferritin_values)

if (!file.exists(paste0("./data/PUBL_finrisk_GlycA_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "FinRisk97", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        bootobs$men[[i]] <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        # Women|Menstr
        bootobs$women_mens[[i]] <- boot(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        # Women|Non-menstr
        bootobs$women_nonmens[[i]] <- boot(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                    var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        
        ## Store
        # Men
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Menstr
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        # Women|Non-menstr
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_mens, means_women_nonmens),
                                     upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                     lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                     Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_finrisk_GlycA_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_finrisk_GlycA_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    # Combine
    means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

if (!file.exists(paste0("./data/PUBL_health2k_GlycA_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    GlycA_trld <- median(fer_glyca[fer_glyca$Cohort == "Health2k", "GlycA"][[1]], na.rm = T)
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    bootobs$men[[i]] <- boot(fer_glyca %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
    # Women|Menstr
    bootobs$women_mens[[i]] <- boot(fer_glyca %>% filter(Group == "Women|Menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
    # Women|Non-menstr
    bootobs$women_nonmens[[i]] <- boot(fer_glyca %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = GlycA, var1_trld = ferritin_values[i], var2_trld = GlycA_trld, var2_over = T)
    ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
    
    ## Store
    # Men
    means_men[i] <- bootobs$men[[i]]$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Menstr
    means_women_mens[i] <- bootobs$women_mens[[i]]$t0
    upper_women_mens[i] <- ci_obj_women_mens$normal[3]
    lower_women_mens[i] <- ci_obj_women_mens$normal[2]
    # Women|Non-menstr
    means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
    upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
    lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))    
    
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_health2k_GlycA_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_health2k_GlycA_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 153), rep("Health2k", 153))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Menstr", "Women|Non-menstr", "Men"))

Proportion analysis plot (GlycA)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") + 
    theme(legend.position = "none")

Ferritin X HbA1C | COLORED

options(scipen = 10000)
ratio4 <- round(nrow(fer_hba1c %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", HbA1C >= 42 & Ferritin >= 15)) * 100 / nrow(fer_hba1c %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_hba1c %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", HbA1C >= 42 & Ferritin >= 15)) * 100 / nrow(fer_hba1c %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_hba1c %>% filter(Group == "Men" & Cohort == "Health2k", HbA1C >= 42 & Ferritin >= 15)) * 100 / nrow(fer_hba1c %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 3), HbA1C = rep(60, 3), lab = c(paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = rep("Health2k", 3), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_hba1c, aes(x = Ferritin, y = HbA1C)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 42, ymax = 100, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group)) +
    theme(legend.position = "none") +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = expression(paste(HbA[1*C], " (mmol/mol)")))

Proportion analysis code (HbA1C)

The reference values for healthy people are between 20 - 42 mmol/mol (https://www.terveyskirjasto.fi/snk03092). Because the relationship with ferritin appears to be positive, we’ll use the upper bound.

iterations <- length(ferritin_values)
HbA1C_trld <- 42

if (!file.exists(paste0("./data/PUBL_health2k_HbA1C_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    bootobs$men[[i]] <- boot(fer_hba1c %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
    # Women|Menstr
    bootobs$women_mens[[i]] <- boot(fer_hba1c %>% filter(Group == "Women|Menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
    # Women|Non-menstr
    bootobs$women_nonmens[[i]] <- boot(fer_hba1c %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                                var1 = Ferritin, var2 = HbA1C, var1_trld = ferritin_values[i], var2_trld = HbA1C_trld, var2_over = T)
    ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
    
    ## Store
    # Men
    means_men[i] <- bootobs$men[[i]]$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Menstr
    means_women_mens[i] <- bootobs$women_mens[[i]]$t0
    upper_women_mens[i] <- ci_obj_women_mens$normal[3]
    lower_women_mens[i] <- ci_obj_women_mens$normal[2]
    # Women|Non-menstr
    means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
    upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
    lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))    
    
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_health2k_HbA1C_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_health2k_HbA1C_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }    
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

means_all <- means_health2k
means_all$Cohort <- c(rep("Health2k", 153))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Menstr", "Women|Non-menstr", "Men"))

Proportion analysis plot (HbA1C)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") + 
    theme(legend.position = "none")

Ferritin X APOB | COLORED

options(scipen = 10000)
ratio1 <- round(nrow(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97")), 2)
ratio3 <- round(nrow(fer_apob %>% filter(Group == "Men" & Cohort == "FinRisk97", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Men" & Cohort == "FinRisk97")), 2)
ratio4 <- round(nrow(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_apob %>% filter(Group == "Men" & Cohort == "Health2k", APOB >= 1.3 & Ferritin >= 15)) * 100 / nrow(fer_apob %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 6), APOB = rep(2.2, 6), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%"), paste0(ratio3, "%"), paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = c(rep("FinRisk97", 3), rep("Health2k", 3)), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_apob, aes(x = Ferritin, y = APOB)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 1.3, ymax = 3, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none") +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "ApoB (g/l)")

Ferritin X APOA1 | COLORED

options(scipen = 10000)
ratio1 <- round(nrow(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97")), 2)
ratio2 <- round(nrow(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97")), 2)
ratio3 <- round(nrow(fer_apoa1 %>% filter(Group == "Men" & Cohort == "FinRisk97", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Men" & Cohort == "FinRisk97")), 2)
ratio4 <- round(nrow(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "Health2k", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "Health2k")), 2)
ratio5 <- round(nrow(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k")), 2)
ratio6 <- round(nrow(fer_apoa1 %>% filter(Group == "Men" & Cohort == "Health2k", APOA1 <= 1.2 & Ferritin >= 15)) * 100 / nrow(fer_apoa1 %>% filter(Group == "Men" & Cohort == "Health2k")), 2)
ann_text <- data.frame(Ferritin = rep(900, 6), APOA1 = rep(0.6, 6), lab = c(paste0(ratio1, "%"), paste0(ratio2, "%"), paste0(ratio3, "%"), paste0(ratio4, "%"), paste0(ratio5, "%"), paste0(ratio6, "%")), Cohort = c(rep("FinRisk97", 3), rep("Health2k", 3)), Group = rep(c("Women|Menstr", "Women|Non-menstr", "Men"), 2))
ggplot(data = fer_apoa1, aes(x = Ferritin, y = APOA1)) + 
    annotate("rect", xmin = 15, xmax = 2000, ymin = 0.4, ymax = 1.2, alpha = .5, fill = "grey") +
    geom_text(data = ann_text, aes(label = lab)) +
    geom_point(aes(color = Group), alpha = 0.1) +
    scale_x_log10() +
    scale_y_log10() +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() + 
    geom_smooth(method = "lm", color = "black", linetype = "dashed", size = 0.5) +
    stat_cor(aes(label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), p.accuracy = 0.001) +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    theme(legend.position = "none") +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "ApoA1 (g/l)")

Proportion analysis code (APOB)

Apolipoprotein B elevates with the risk of cardiovascular diseases and it correlates positively with ferritin. The reference values for healthy men are 0.6 - 1.5 (of which we’ll use 1.5) and for healthy women they are 0.6 - 1.3 (we’ll use 1.3). (https://huslab.fi/ohjekirja/20705.html)

iterations <- length(ferritin_values)
APOB_trld_men <- 1.5 # for males [0.6, 1.5]
APOB_trld_women <- 1.3 # for females [0.6, 1.3]

if (!file.exists(paste0("./data/PUBL_finrisk_APOB_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        bootobs$men[[i]] <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        # Women|Menstr
        bootobs$women_mens[[i]] <- boot(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        # Women|Non-menstr
        bootobs$women_nonmens[[i]] <- boot(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        
        ## Store
        # Men
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Menstr
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        # Women|Non-menstr
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_mens, means_women_nonmens),
                                     upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                     lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                     Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_finrisk_APOB_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_finrisk_APOB_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    # Combine
    means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

if (!file.exists(paste0("./data/PUBL_health2k_APOB_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    bootobs$men[[i]] <- boot(fer_apob %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_men, var2_over = T)
    ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
    # Women|Menstr
    bootobs$women_mens[[i]] <- boot(fer_apob %>% filter(Group == "Women|Menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
    # Women|Non-menstr
    bootobs$women_nonmens[[i]] <- boot(fer_apob %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOB, var1_trld = ferritin_values[i], var2_trld = APOB_trld_women, var2_over = T)
    ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
    
    ## Store
    # Men
    means_men[i] <- bootobs$men[[i]]$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Menstr
    means_women_mens[i] <- bootobs$women_mens[[i]]$t0
    upper_women_mens[i] <- ci_obj_women_mens$normal[3]
    lower_women_mens[i] <- ci_obj_women_mens$normal[2]
    # Women|Non-menstr
    means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
    upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
    lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))    
    
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_health2k_APOB_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_health2k_APOB_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 153), rep("Health2k", 153))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Menstr", "Women|Non-menstr", "Men"))

Proportion analysis plot (APOB)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") + 
    theme(legend.position = "none")

Proportion analysis code (APOA1)

Apolipoprotein A1 decreases with the risk of atherosclerosis (peripheral, coronary, brain) and appears to correlate negatively with ferritin. Thus, we will use the lower bounds of the reference intervals for healthy individuals here. These are 1.1 - 2.0 for men and 1.2 - 2.3 for women. (https://huslab.fi/ohjekirja/20705.html)

iterations <- length(ferritin_values)
APOA1_trld_men <- 1.1 # for males [1.1, 2.0]
APOA1_trld_women <- 1.2 # for females [1.2, 2.3]

if (!file.exists(paste0("./data/PUBL_finrisk_APOA1_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))
    
    for (i in 1:iterations) {
        
        #############
        #### FinRisk97
        #############
        
        ## Compute
        # Men
        bootobs$men[[i]] <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                             var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        # Women|Menstr
        bootobs$women_mens[[i]] <- boot(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        # Women|Non-menstr
        bootobs$women_nonmens[[i]] <- boot(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "FinRisk97"), statistic = get_ratio_boot, R = boot_n, 
                                   var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        
        ## Store
        # Men
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        # Women|Menstr
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        # Women|Non-menstr
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
        
        # Combine
        means_finrisk <- data.frame(Ferritin = rep(ferritin_values, 3),
                                     means = c(means_men, means_women_mens, means_women_nonmens),
                                     upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                     lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                     Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_finrisk_APOA1_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_finrisk_APOA1_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

if (!file.exists(paste0("./data/PUBL_health2k_APOA1_", boot_n, ".rds"))) { # run bootstrap only if needed
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    bootobs <- setNames(vector("list", length = 3), c("men", "women_mens", "women_nonmens"))

    for (i in 1:iterations) {
    
    #############
    #### Health2000
    #############
    
    ## Compute
    # Men
    bootobs$men[[i]] <- boot(fer_apoa1 %>% filter(Group == "Men" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                         var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_men, var2_over = F)
    ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
    # Women|Menstr
    bootobs$women_mens[[i]] <- boot(fer_apoa1 %>% filter(Group == "Women|Menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
    # Women|Non-menstr
    bootobs$women_nonmens[[i]] <- boot(fer_apoa1 %>% filter(Group == "Women|Non-menstr" & Cohort == "Health2k"), statistic = get_ratio_boot, R = boot_n, 
                               var1 = Ferritin, var2 = APOA1, var1_trld = ferritin_values[i], var2_trld = APOA1_trld_women, var2_over = F)
    ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
    
    ## Store
    # Men
    means_men[i] <- bootobs$men[[i]]$t0
    upper_men[i] <- ci_obj_men$normal[3]
    lower_men[i] <- ci_obj_men$normal[2]
    # Women|Menstr
    means_women_mens[i] <- bootobs$women_mens[[i]]$t0
    upper_women_mens[i] <- ci_obj_women_mens$normal[3]
    lower_women_mens[i] <- ci_obj_women_mens$normal[2]
    # Women|Non-menstr
    means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
    upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
    lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))    
    
    
    }
    
    # Save
    saveRDS(bootobs, paste0("./data/PUBL_health2k_APOA1_", boot_n, ".rds"))
} else {
    bootobs <- readRDS(paste0("./data/PUBL_health2k_APOA1_", boot_n, ".rds"))
    
    ## Preallocate
    # Men
    means_men <- 1:iterations
    upper_men <- 1:iterations
    lower_men <- 1:iterations
    # Women|Menstr
    means_women_mens <- 1:iterations
    upper_women_mens <- 1:iterations
    lower_women_mens <- 1:iterations
    # Women|Non-menstr
    means_women_nonmens <- 1:iterations
    upper_women_nonmens <- 1:iterations
    lower_women_nonmens <- 1:iterations
    
    for (i in 1:iterations) {
        ci_obj_men <- boot.ci(bootobs$men[[i]], type = "norm")
        ci_obj_women_mens <- boot.ci(bootobs$women_mens[[i]], type = "norm")
        ci_obj_women_nonmens <- boot.ci(bootobs$women_nonmens[[i]], type = "norm")
        # Store
        means_men[i] <- bootobs$men[[i]]$t0
        upper_men[i] <- ci_obj_men$normal[3]
        lower_men[i] <- ci_obj_men$normal[2]
        means_women_mens[i] <- bootobs$women_mens[[i]]$t0
        upper_women_mens[i] <- ci_obj_women_mens$normal[3]
        lower_women_mens[i] <- ci_obj_women_mens$normal[2]
        means_women_nonmens[i] <- bootobs$women_nonmens[[i]]$t0
        upper_women_nonmens[i] <- ci_obj_women_nonmens$normal[3]
        lower_women_nonmens[i] <- ci_obj_women_nonmens$normal[2]
    }    
    
    # Combine
    means_health2k <- data.frame(Ferritin = rep(ferritin_values, 3),
                                 means = c(means_men, means_women_mens, means_women_nonmens),
                                 upper = c(upper_men, upper_women_mens, upper_women_nonmens),
                                 lower = c(lower_men, lower_women_mens, lower_women_nonmens),
                                 Gender = c(rep("Men", iterations), rep("Women|Menstr", iterations), rep("Women|Non-menstr", iterations)))
    }

means_all <- rbind(means_finrisk, means_health2k)
means_all$Cohort <- c(rep("FinRisk97", 153), rep("Health2k", 153))
means_all$Group <- factor(means_all$Gender, levels = c("Women|Menstr", "Women|Non-menstr", "Men"))

Proportion analysis plot (APOA1)

ggplot(data = means_all, aes(x = Ferritin, y = means)) +
    geom_ribbon(aes(ymin = lower, ymax = upper, fill = Group), alpha = .3) +
    geom_line(aes(linetype = Group)) +
    scale_color_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                       limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    scale_fill_manual(values = c( "#00BFFF",  "#de2d26", "#ff85a2" ),
                      limits = c( "Men",  "Women|Non-menstr", "Women|Menstr" )) +
    theme_minimal() +
    facet_grid(rows = vars(Group), cols = vars(Cohort)) +
    labs(x = expression(paste("Ferritin (", mu, "g/l)")), 
         y = "Change (pp)") + 
    theme(legend.position = "none")